R code for ‘How changes in spring and autumn phenology translate into growth - experimental evidence of asymmetric effects’
1. Load packages
get packages
require(ggplot2)
require(dplyr)
require(patchwork)
require(lubridate)
require(data.table)
require(git2r)
require(stringr)
require(broom)
require(viridis)
require(car)2. Load data
data tables
# define paths
data.dir = "/Users/consti/Desktop/PhD/Publication_material/10_PhenogrowthEx3/Rdata/Experiment"
output.dir = "/Users/consti/Desktop/PhD/Publication_material/10_PhenogrowthEx3/Routput"
#load data
clim.data = data.frame(fread(paste(data.dir,"climate_data.csv",sep="/"), header=T, sep=",")) # Climate
volume.data = data.frame(fread(paste(data.dir,"volume.csv",sep="/"), header=T, sep=",")) # Stem volume
pheno.data = data.frame(fread(paste(data.dir,"phenology.csv",sep="/"))) # Phenology
biomass.data = data.frame(fread(paste(data.dir,"biomass.csv",sep="/"))) # Biomassggplot themes
#define plot themes
plotTheme = theme(legend.position = "right",
legend.background = element_rect(fill=NA, size=0.5, linetype="solid"),
legend.text = element_text(color="black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(colour = NA, size=1, fill="grey97"),
axis.line = element_line(color = "black"),
strip.background = element_rect(fill=NA),
strip.text = element_text(colour = 'black'))
plotTheme2 = theme(legend.position = "none",
legend.background = element_rect(fill=NA, size=0.5, linetype="solid"),
legend.text = element_text(color="black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(colour = NA, size=1, fill="grey97"),
axis.text = element_blank(),
axis.ticks = element_blank(),
strip.background = element_rect(fill=NA),
strip.text = element_text(colour = 'black'))
plotTheme3 = theme(legend.position = "none",
legend.background = element_rect(fill=NA, size=0.5, linetype="solid"),
legend.text = element_text(color="black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(colour = NA, size=1, fill="grey97"),
axis.line = element_line(color = "black"),
strip.background = element_rect(fill=NA),
strip.text = element_text(colour = 'black'))
plotTheme4 = theme(legend.position = "right",
legend.background = element_rect(fill=NA, size=0.5, linetype="solid"),
legend.text = element_text(color="black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(colour = NA, size=1, fill="grey97"),
axis.text.x = element_text(angle = 45, hjust = 1),
axis.line.y = element_line(color = "black"),
strip.background = element_rect(fill=NA),
strip.text = element_text(colour = 'black'))3. Data preparation
climate data
data management
#Convert dates to DOY
clim.data$date = dmy(clim.data$date)
#get median deviations from control + interquartile ranges
clim.data1 = data.frame(clim.data %>%
filter(!is.na(inside_tmp)) %>% #delete summer period
mutate(diffOutIn = inside_tmp-outside_tmp)) %>% #add deviation column
group_by(treatment) %>%
summarize(medianDiff = median(diffOutIn),
IQrange = IQR(diffOutIn, na.rm = FALSE, type = 7))
#get new dataframe for deviation plot
clim.data2 = data.frame(clim.data %>%
filter(!is.na(inside_tmp)) %>% #delete NAs
mutate(diffOutIn = inside_tmp-outside_tmp)) #get deviation from control
#reshape initial dataframe
clim.data3 = melt(clim.data, id.vars = c("date","treatment"), measure.vars = c("outside_tmp","inside_tmp"))
#create treatment column
clim.data3$variable2 = ifelse(clim.data3$variable=="outside_tmp", "control",
ifelse(clim.data$treatment=="autumn", "autumn", "spring"))temperature table
data.frame(clim.data1)## treatment medianDiff IQrange
## 1 autumn 4.6715 1.94125
## 2 spring 4.7170 2.23475stem volume
stem volume calculation
#order treatments
Treatment = c("spring", "autumn", "autumn-spring", "control")
volume.data$Treatment = factor(volume.data$Treatment, levels=Treatment)
#calculate volumes
volume.data$Volume = 1/3 * pi * volume.data$Height *
((volume.data$StemDiameter/2)^2 +
(volume.data$StemDiameter/2)*(volume.data$StemDiameterTop/2) +
volume.data$StemDiameterTop^2)
#sum stem and twig volume measurements within individuals for each year
VariablesVector = c("Volume")
volume.data = data.frame(volume.data %>%
group_by(Species,ID,Treatment,Year) %>%
summarize_at(VariablesVector, sum, na.rm = TRUE))sample sizes
table(volume.data$Treatment, volume.data$Species)/4##
## Fagus sylvatica Lonicera xylosteum Quercus robur
## spring 9 10 8
## autumn 9 8 10
## autumn-spring 9 9 10
## control 10 9 9phenology
#get RMFs
biomass.data$RSR = biomass.data$Root/biomass.data$Shoot
#delete species epithet in all tables
DeleteEpithet <- function(df) {
df$Species = gsub(' [A-z ]*', '' , df$Species)
return(df)
}
dfnames = list(pheno.data,biomass.data,volume.data)
dfs = lapply(dfnames, DeleteEpithet)
names(dfs) = c("pheno.data","biomass.data","volume.data")
#Ordering of factors
Treatments = c("control", "autumn", "autumn-spring", "spring") #order treatments
pheno.data$Treatment = factor(pheno.data$Treatment, levels=Treatments, ordered=T)
Species = c("Quercus", "Fagus", "Lonicera") #order species
pheno.data$Species = factor(pheno.data$Species, levels=Species, ordered=T)
#reshape phenology table to long format
pheno.data = melt(dfs$pheno.data, id.vars = c("Species","ID","Treatment"))
#Convert dates to DOY
pheno.data$value = dmy(pheno.data$value)
pheno.data$DOY = yday(pheno.data$value) # Jan 1 = day 1
#create year and phenophase columns
pheno.data$Year = str_sub(pheno.data$variable, str_length(pheno.data$variable)-3, str_length(pheno.data$variable))
pheno.data$Pheno = str_sub(pheno.data$variable, 1, str_length(pheno.data$variable)-5)
pheno.data = select(pheno.data, -c(variable, value))
pheno.data = na.omit(pheno.data)
#from long (many rows) to short (multiple columns) format
pheno.data = dcast(pheno.data, Species + ID + Treatment + Year ~ Pheno, value.var = "DOY")
#get growing season length information (leaf-out to senescence)
pheno.data$GSL = pheno.data$Senescence - pheno.data$Leaf_out
#constrain late leafers
pheno.data$Leaf_out = ifelse(pheno.data$Leaf_out>140,140,pheno.data$Leaf_out)biomass
#Melt data to a long format
biomass.data = melt(dfs$biomass.data, id = c("Species", "ID", "Treatment"),
variable.name="organ", value.name="biomass")
#Ordering of factors
Treatments = c("control", "autumn", "autumn-spring", "spring") #order treatments
biomass.data$Treatment = factor(biomass.data$Treatment, levels=Treatments, ordered=T)
Species = c("Quercus", "Fagus", "Lonicera") #order species
biomass.data$Species = factor(biomass.data$Species, levels=Species, ordered=T)
variable = c("Total", "Shoot", "Root", "RSR") #order organ
biomass.data$organ = factor(biomass.data$organ, levels=variable, ordered=T)4. Analysis
phenology
treatment effect on phenology
#################################
# Summarize annual phenology data
#################################
#get average phenology anomalies per individual
###############################################
#reshape phenology columns to long format
pheno.data = melt(pheno.data, id.vars = c("Species","ID","Treatment","Year"),
variable.name="phenology", value.name="day")
#get means per species and year
phenoMean = pheno.data %>%
group_by(Species,phenology,Year) %>%
summarize(Mean = mean(day, na.rm=TRUE),
SD = sd(day, na.rm=TRUE))
#merge data
pheno.data <-merge(pheno.data, phenoMean, all.x=T, by=c("Species","phenology","Year"))
#get change in phenology relative to mean
pheno.data = pheno.data %>%
mutate(phenologyAnomaly = ifelse(phenology=="Leaf_out", Mean-day, day-Mean)) %>% #add column
dplyr::select(!c(Mean)) #delete columns
#from long (many rows) to short (multiple columns) format
pheno.data3 = dcast(pheno.data, Species + ID + Treatment + Year ~ phenology, value.var = c("phenologyAnomaly"))
#summarize years
VariablesVector = c("Leaf_out","Senescence","GSL")
pheno.data3 = data.frame(pheno.data3 %>%
filter(!Year %in% c("2017")) %>% #delete rows based on condition
group_by(Species,ID,Treatment) %>%
summarize_at(VariablesVector, mean, na.rm = TRUE))
# Add biomass information
#########################
#merge phenology and biomass tables
pheno.data3 = merge(pheno.data3, dfs$biomass.data, all.x=T, by=c("Species","ID","Treatment"))
#reshape phenology columns to long format
pheno.data3 = melt(pheno.data3, id.vars = c("Species","ID","Treatment","Total","Root","Shoot","RSR"),
variable.name="phenology", value.name="phenologyAnomaly")
#get average phenology dates per individual
###########################################
#from long (many rows) to short (multiple columns) format
pheno.data1 = dcast(pheno.data, Species + ID + Treatment + Year ~ phenology, value.var = c("day"))
#summarize years
VariablesVector = c("Leaf_out","Senescence","GSL")
pheno.data1 = data.frame(pheno.data1 %>%
filter(!Year %in% c("2017")) %>% #delete rows based on condition
group_by(Species,ID,Treatment) %>%
summarize_at(VariablesVector, mean, na.rm = TRUE))
#merge phenology and biomass tables
pheno.data1 = merge(pheno.data1, dfs$biomass, all.x=T, by=c("Species","ID","Treatment"))
#reshape phenology columns to long format
pheno.data1 = melt(pheno.data1, id.vars = c("Species","ID","Treatment","Total","Root","Shoot","RSR"),
variable.name="phenology", value.name="day")
#########################################################################################
#########################################################################################
###########################
# Multivariate linear model
###########################
#Create spring and autumn treatments
####################################
pheno.data4 = pheno.data3
pheno.data4$spring = ifelse(pheno.data4$Treatment %in% c("control","autumn"), "no","warming")
pheno.data4$autumn = ifelse(pheno.data4$Treatment %in% c("control","spring"), "no","warming")
#Multivariate linear model
##########################
resultsLM.pheno = pheno.data4 %>%
group_by(Species,phenology) %>%
do({model = lm(phenologyAnomaly ~ spring*autumn, data=.) # create your model
data.frame(tidy(model))}) %>% # get coefficient info
filter(!term %in% c("(Intercept)")) %>% # delete rows
dplyr::select(Species, phenology, term, p.value, estimate, std.error) %>% #delete columns
mutate_if(is.numeric, round, 2) %>%
mutate(sig = ifelse(p.value<0.05, "**", ifelse(p.value<=0.1, "*", ""))) #add significance
#Ordering of factors
term = c("springwarming", "autumnwarming", "springwarming:autumnwarming") #order treatments
resultsLM.pheno$term = factor(resultsLM.pheno$term, levels=term, ordered=T)
Species = c("Quercus", "Fagus", "Lonicera") #order species
resultsLM.pheno$Species = factor(resultsLM.pheno$Species, levels=Species, ordered=T)
#order statistics table to match dataframe
resultsLM.pheno = resultsLM.pheno[with(resultsLM.pheno, order(resultsLM.pheno$Species, resultsLM.pheno$phenology)), ]
data.frame(resultsLM.pheno[,-c(7)])## Species phenology term p.value estimate std.error
## 1 Quercus Leaf_out springwarming 0.00 16.27 1.16
## 2 Quercus Leaf_out autumnwarming 0.06 -2.17 1.09
## 3 Quercus Leaf_out springwarming:autumnwarming 0.79 0.43 1.57
## 4 Quercus Senescence springwarming 0.02 -10.77 4.56
## 5 Quercus Senescence autumnwarming 0.00 13.67 4.31
## 6 Quercus Senescence springwarming:autumnwarming 0.39 -5.43 6.19
## 7 Quercus GSL springwarming 0.27 5.50 4.88
## 8 Quercus GSL autumnwarming 0.02 11.50 4.61
## 9 Quercus GSL springwarming:autumnwarming 0.46 -5.00 6.63
## 10 Fagus Leaf_out springwarming 0.00 13.72 1.66
## 11 Fagus Leaf_out autumnwarming 0.00 -6.73 1.66
## 12 Fagus Leaf_out springwarming:autumnwarming 0.71 0.89 2.37
## 13 Fagus Senescence springwarming 0.12 -3.28 2.05
## 14 Fagus Senescence autumnwarming 0.00 23.11 2.05
## 15 Fagus Senescence springwarming:autumnwarming 0.74 -1.00 2.94
## 16 Fagus GSL springwarming 0.00 10.44 2.21
## 17 Fagus GSL autumnwarming 0.00 16.38 2.21
## 18 Fagus GSL springwarming:autumnwarming 0.97 -0.11 3.16
## 19 Lonicera Leaf_out springwarming 0.00 21.88 1.95
## 20 Lonicera Leaf_out autumnwarming 0.09 -3.54 2.06
## 21 Lonicera Leaf_out springwarming:autumnwarming 0.00 -8.84 2.83
## 22 Lonicera Senescence springwarming 0.04 -4.75 2.19
## 23 Lonicera Senescence autumnwarming 0.00 7.69 2.32
## 24 Lonicera Senescence springwarming:autumnwarming 0.08 5.67 3.19
## 25 Lonicera GSL springwarming 0.00 17.13 2.77
## 26 Lonicera GSL autumnwarming 0.17 4.15 2.93
## 27 Lonicera GSL springwarming:autumnwarming 0.44 -3.17 4.04######################################################
# Check for deviation from control (one-sample T-test)
######################################################
#get means of control group
pheno.data2 = pheno.data1 %>%
filter(Treatment %in% c("control")) %>%
group_by(Species,phenology) %>%
summarize(Mean = mean(day, na.rm=TRUE))
#merge data
pheno.data2 <-merge(pheno.data1[!pheno.data1$Treatment=="control",],
pheno.data2, all.x=T, by=c("Species","phenology"))
#get change in phenology relative to control
pheno.data2$phenologyChange = pheno.data2$day - pheno.data2$Mean
#Ordering of factors
Phenology = c("Leaf_out", "Senescence", "GSL") #order treatments
pheno.data2$phenology = factor(pheno.data2$phenology, levels=Phenology, ordered=T)
Treatments = c("spring", "autumn", "autumn-spring") #order treatments
pheno.data2$Treatment = factor(pheno.data2$Treatment, levels=Treatments, ordered=T)
Species = c("Quercus", "Fagus", "Lonicera") #order species
pheno.data2$Species = factor(pheno.data2$Species, levels=Species, ordered=T)
#T-test
resultsTT = pheno.data2 %>%
group_by(Species,phenology,Treatment) %>%
do({model = t.test(.$phenologyChange) # create your model
data.frame(tidy(model),
tidy(shapiro.test(x = .$phenologyChange))[1,2])}) %>%
dplyr::select(Species, phenology, Treatment, p.value, estimate, p.value.1) %>% #delete columns
rename(shapiroTest=p.value.1) %>%
mutate(significance = ifelse(p.value<0.05, "**", ifelse(p.value<=0.1, "*", ""))) %>% #add significance
mutate_if(is.numeric, round, 2)
#order statistics table to match dataframe
resultsTT = resultsTT[with(resultsTT, order(resultsTT$Species, resultsTT$phenology)), ]
data.frame(resultsTT[,-c(7)])## Species phenology Treatment p.value estimate shapiroTest
## 1 Quercus Leaf_out spring 0.00 -16.27 0.41
## 2 Quercus Leaf_out autumn 0.01 2.17 0.65
## 3 Quercus Leaf_out autumn-spring 0.00 -14.53 0.24
## 4 Quercus Senescence spring 0.00 -10.77 0.36
## 5 Quercus Senescence autumn 0.00 13.67 0.73
## 6 Quercus Senescence autumn-spring 0.52 -2.53 0.37
## 7 Quercus GSL spring 0.03 5.50 0.33
## 8 Quercus GSL autumn 0.01 11.50 0.69
## 9 Quercus GSL autumn-spring 0.01 12.00 0.38
## 10 Fagus Leaf_out spring 0.00 -13.72 0.41
## 11 Fagus Leaf_out autumn 0.00 6.73 0.17
## 12 Fagus Leaf_out autumn-spring 0.00 -7.88 0.77
## 13 Fagus Senescence spring 0.08 -3.28 0.49
## 14 Fagus Senescence autumn 0.00 23.11 0.02
## 15 Fagus Senescence autumn-spring 0.00 18.83 0.22
## 16 Fagus GSL spring 0.00 10.44 0.60
## 17 Fagus GSL autumn 0.00 16.38 0.31
## 18 Fagus GSL autumn-spring 0.00 26.72 0.42
## 19 Lonicera Leaf_out spring 0.00 -21.88 0.08
## 20 Lonicera Leaf_out autumn 0.01 3.54 0.36
## 21 Lonicera Leaf_out autumn-spring 0.00 -9.50 0.98
## 22 Lonicera Senescence spring 0.02 -4.75 0.48
## 23 Lonicera Senescence autumn 0.00 7.69 0.02
## 24 Lonicera Senescence autumn-spring 0.00 8.61 0.01
## 25 Lonicera GSL spring 0.00 17.13 0.44
## 26 Lonicera GSL autumn 0.00 4.15 0.36
## 27 Lonicera GSL autumn-spring 0.00 18.11 0.98phenology effect on total biomass
##########################
# Univariate linear models
##########################
#reshape biomass columns to long format
pheno.data3 = melt(pheno.data3, id.vars = c("Species","ID","Treatment","phenology","phenologyAnomaly"),
variable.name="organ", value.name="biomass")
#Transform biomass to percentage anomaly
########################################
#get means per species
pheno.data5 = pheno.data3 %>%
group_by(Species,organ) %>%
summarize(Mean = mean(biomass, na.rm=TRUE))
#merge data
pheno.data3 <-merge(pheno.data3, pheno.data5, all.x=T, by=c("Species","organ"))
#get precentage change in biomass
pheno.data3$relativeBiomass = (pheno.data3$biomass / pheno.data3$Mean -1) * 100
#order species
Species = c("Quercus", "Fagus", "Lonicera")
pheno.data3$Species = factor(pheno.data3$Species, levels=Species, ordered=T)
#Extract linear model info
resultsLM.pheno2 = pheno.data3 %>%
group_by(Species,phenology,organ) %>%
do({model = lm(relativeBiomass ~ phenologyAnomaly, data=.) # create model
data.frame(tidy(model), # get coefficient info
glance(model),
tidy(shapiro.test(x = residuals(object = model)))[1,2]
)}) %>% # get model info
filter(term != "(Intercept)") %>% # delete intercept info
dplyr::select(Species, phenology, organ, estimate, std.error, p.value, r.squared, p.value.2) %>% # delete columns
rename(shapiro=p.value.2) %>% #rename columns
mutate(sig = ifelse(p.value<0.05, "**", ifelse(p.value<=0.1, "*", ""))) %>% # add column
mutate_if(is.numeric, round, 2)
#Order phenophases
Organs = c("Total","Shoot","Root","RSR") #order treatments
resultsLM.pheno2$organ = factor(resultsLM.pheno2$organ, levels=Organs, ordered=T)
#order statistics table to match dataframe
resultsLM.pheno2 = resultsLM.pheno2[with(resultsLM.pheno2,
order(resultsLM.pheno2$Species,
resultsLM.pheno2$organ,
resultsLM.pheno2$phenology)), ]
data.frame(resultsLM.pheno2)## Species phenology organ estimate std.error p.value r.squared shapiro sig
## 1 Quercus Leaf_out Total 2.46 0.77 0.00 0.22 0.39 **
## 2 Quercus Senescence Total -2.12 0.49 0.00 0.35 0.95 **
## 3 Quercus GSL Total -1.30 0.67 0.06 0.10 0.35 *
## 4 Quercus Leaf_out Shoot 1.79 0.74 0.02 0.14 0.23 **
## 5 Quercus Senescence Shoot -1.98 0.44 0.00 0.37 0.61 **
## 6 Quercus GSL Shoot -1.53 0.59 0.01 0.16 0.22 **
## 7 Quercus Leaf_out Root 3.03 0.89 0.00 0.25 0.62 **
## 8 Quercus Senescence Root -2.25 0.60 0.00 0.29 0.75 **
## 9 Quercus GSL Root -1.10 0.81 0.18 0.05 0.02
## 10 Quercus Leaf_out RSR 0.98 0.54 0.08 0.09 0.83 *
## 11 Quercus Senescence RSR -0.18 0.39 0.64 0.01 0.16
## 12 Quercus GSL RSR 0.38 0.45 0.40 0.02 0.36
## 13 Fagus Leaf_out Total 0.78 0.51 0.14 0.06 0.28
## 14 Fagus Senescence Total -0.16 0.36 0.66 0.01 0.60
## 15 Fagus GSL Total 0.27 0.41 0.52 0.01 0.44
## 16 Fagus Leaf_out Shoot 0.42 0.52 0.43 0.02 0.31
## 17 Fagus Senescence Shoot 0.09 0.36 0.80 0.00 0.34
## 18 Fagus GSL Shoot 0.37 0.40 0.36 0.02 0.50
## 19 Fagus Leaf_out Root 1.24 0.60 0.05 0.11 0.18 **
## 20 Fagus Senescence Root -0.49 0.43 0.27 0.03 0.57
## 21 Fagus GSL Root 0.13 0.50 0.80 0.00 0.07
## 22 Fagus Leaf_out RSR 0.81 0.42 0.06 0.10 0.28 *
## 23 Fagus Senescence RSR -0.50 0.29 0.10 0.08 0.82
## 24 Fagus GSL RSR -0.14 0.34 0.68 0.00 0.94
## 25 Lonicera Leaf_out Total -0.16 0.20 0.45 0.02 0.39
## 26 Lonicera Senescence Total 0.37 0.30 0.22 0.04 0.44
## 27 Lonicera GSL Total 0.01 0.22 0.96 0.00 0.44
## 28 Lonicera Leaf_out Shoot -0.13 0.22 0.58 0.01 0.75
## 29 Lonicera Senescence Shoot 0.07 0.34 0.84 0.00 0.90
## 30 Lonicera GSL Shoot -0.12 0.25 0.64 0.01 0.79
## 31 Lonicera Leaf_out Root -0.21 0.30 0.50 0.01 0.50
## 32 Lonicera Senescence Root 0.91 0.43 0.04 0.12 0.09 **
## 33 Lonicera GSL Root 0.24 0.33 0.48 0.02 0.18
## 34 Lonicera Leaf_out RSR -0.09 0.33 0.79 0.00 0.27
## 35 Lonicera Senescence RSR 0.87 0.47 0.07 0.09 0.43 *
## 36 Lonicera GSL RSR 0.37 0.36 0.32 0.03 0.26#Visually inspect model assumptions
data.assumptions = pheno.data3
data.assumptions$category = paste(data.assumptions$Species,
data.assumptions$organ,
data.assumptions$phenology, sep="_") #create identifier column
category.list = as.factor(unique(data.assumptions$category)) #create category vector
par(mfrow=c(2,2)) #set plot layout
for (category in category.list){ #loop over categories
tab.subset=data.assumptions[data.assumptions$category==category, ] #create table subset
plot(lm(relativeBiomass ~ phenologyAnomaly, data=tab.subset), main=category)
}
biomass
total and per treatment biomass distributions
#Species-level plot
ggplot(data = biomass.data, mapping = aes(biomass, fill=organ)) +
geom_histogram(bins=10, position="identity", color="black") +
xlab("Biomass (g) [panels 1-3]
Biomass ratio [panel 4]") +
ylab("Count") +
scale_fill_viridis(option="E",discrete=TRUE) +
facet_grid(Species~organ, scales="free") +
plotTheme4
#Treatment-level plot
ggplot(data = biomass.data[biomass.data$organ %in% c("Shoot","Root"),], mapping = aes(biomass, fill=organ)) +
geom_histogram(bins=10, color="black") +
xlab("Biomass (g)") +
ylab("Count") +
scale_fill_viridis(option="C",discrete=TRUE) +
facet_grid(organ+Species~Treatment, scales="free") +
plotTheme4
treatment analysis
#######
# ANOVA
#######
#Extract model info and check ANOVA assumptions
resultsLM = biomass.data %>%
group_by(Species,organ) %>%
do({model = aov(biomass ~ Treatment, data=.) # create your model
data.frame(tidy(model)[1,], # get coefficient info
glance(model)[1,], # get model info
leveneTest=tidy(leveneTest(model))[1,4], # check Homogeneity of variances
tidy(shapiro.test(x = residuals(object = model)))[1,2] # check for Normality
)}) %>%
#select(Species, organ, p.value, r.squared, adj.r.squared, p.value.2, p.value.3) %>% #delete columns
rename(leveneTest=p.value.2, shapiroTest=p.value.3) %>% #rename columns
mutate(sig = ifelse(p.value<0.05, "**", ifelse(p.value<0.1, "*", "")),
Treatment="control") #add significance and dummy column for plotting
data.frame(resultsLM)## Species organ term df sumsq meansq statistic p.value
## 1 Quercus Total Treatment 3 8.972304e+04 2.990768e+04 6.7781574 0.001097137
## 2 Quercus Shoot Treatment 3 1.352895e+04 4.509648e+03 5.2936506 0.004317898
## 3 Quercus Root Treatment 3 3.498875e+04 1.166292e+04 6.5838558 0.001304986
## 4 Quercus RSR Treatment 3 4.736796e-01 1.578932e-01 1.4440644 0.247684805
## 5 Fagus Total Treatment 3 1.121741e+04 3.739135e+03 1.0439249 0.386050827
## 6 Fagus Shoot Treatment 3 8.885922e+02 2.961974e+02 0.2494104 0.861165544
## 7 Fagus Root Treatment 3 6.091089e+03 2.030363e+03 2.1809050 0.108879735
## 8 Fagus RSR Treatment 3 1.603263e-01 5.344209e-02 1.8326934 0.160419578
## 9 Lonicera Total Treatment 3 3.100722e+02 1.033574e+02 0.2485173 0.861777929
## 10 Lonicera Shoot Treatment 3 1.906361e+02 6.354537e+01 0.3070789 0.820067194
## 11 Lonicera Root Treatment 3 3.923083e+02 1.307694e+02 1.1726019 0.335504179
## 12 Lonicera RSR Treatment 3 7.252300e-02 2.417433e-02 1.7831863 0.170190852
## r.squared adj.r.squared sigma statistic.1 p.value.1 df.1 logLik
## 1 0.38126321 0.325014412 66.4256090 6.7781574 0.001097137 4 -205.639193
## 2 0.32489040 0.263516795 29.1872856 5.2936506 0.004317898 4 -175.212262
## 3 0.37442617 0.317555822 42.0884980 6.5838558 0.001304986 4 -188.755791
## 4 0.11604443 0.035684834 0.3306652 1.4440644 0.247684805 4 -9.438122
## 5 0.08667647 0.003647058 59.8481834 1.0439249 0.386050827 4 -201.781144
## 6 0.02217097 -0.066722576 34.4614360 0.2494104 0.861165544 4 -181.358248
## 7 0.16545943 0.089592102 30.5118457 2.1809050 0.108879735 4 -176.854389
## 8 0.14281440 0.064888438 0.1707642 1.8326934 0.160419578 4 15.012317
## 9 0.02276804 -0.068847458 20.3935328 0.2485173 0.861777929 4 -157.509534
## 10 0.02798305 -0.063143535 14.3852349 0.3070789 0.820067194 4 -144.944976
## 11 0.09904344 0.014578768 10.5603384 1.1726019 0.335504179 4 -133.817484
## 12 0.14322951 0.062907271 0.1164338 1.7831863 0.170190852 4 28.453889
## AIC BIC deviance df.residual leveneTest shapiroTest sig
## 1 421.27839 429.33297 1.456079e+05 33 0.40003230 0.87098652 **
## 2 360.42452 368.47911 2.811262e+04 33 0.98423338 0.99514138 **
## 3 387.51158 395.56617 5.845758e+04 33 0.04824924 0.99099527 **
## 4 28.87624 36.93083 3.608202e+00 33 0.10392883 0.66060872
## 5 413.56229 421.61688 1.181996e+05 33 0.23193665 0.08508579
## 6 372.71650 380.77108 3.919049e+04 33 0.14377313 0.16415740
## 7 363.70878 371.76337 3.072210e+04 33 0.64350555 0.44789637
## 8 -20.02463 -11.97004 9.622935e-01 33 0.49278732 0.52158902
## 9 325.01907 332.93666 1.330868e+04 32 0.92356152 0.53347485
## 10 299.88995 307.80755 6.621919e+03 32 0.70338321 0.69788510
## 11 277.63497 285.55256 3.568664e+03 32 0.77713818 0.48212633
## 12 -46.90778 -38.99018 4.338182e-01 32 0.26694805 0.66541472
## Treatment
## 1 control
## 2 control
## 3 control
## 4 control
## 5 control
## 6 control
## 7 control
## 8 control
## 9 control
## 10 control
## 11 control
## 12 control# Plots to check ANOVA assumptions
biomass.data$category = paste(biomass.data$Species,biomass.data$organ, sep="_") #create identifier column
category.list = as.factor(unique(biomass.data$category)) #create category vector
par(mfrow=c(2,2)) #set plot layout
for (category in category.list){ #loop over categories
tab.subset=biomass.data[biomass.data$category==category, ] #create table subset
plot(aov(biomass ~ Treatment, data=tab.subset), 1, main=category) # 1. Homogeneity of variances
plot(aov(biomass ~ Treatment, data=tab.subset), 2, main=category) # 2. Normality
}
#####
#Plot
#####
Boxp <- ggplot(data = biomass.data, mapping = aes(x = Treatment, y = biomass, fill=Treatment)) +
geom_boxplot(position=position_dodge(0.8),coef=1e30) +
geom_text(data = resultsLM,
mapping = aes(x = Inf, y = Inf, hjust = 2, vjust = 2,
label = paste("R2 = ", round(r.squared,2), sig, sep="")),
size=3.5, color="black")+
xlab("Treatment") +
ylab("Biomass (g)") +
scale_fill_viridis(option="E",discrete=TRUE) +
facet_grid(organ~Species, scales="free") +
plotTheme4
Boxp
Relative biomass analysis
###########################
# Multivariate linear model
###########################
#Create spring and autumn treatments
####################################
biomass.data$spring = ifelse(biomass.data$Treatment %in% c("control","autumn"), "no","warming")
biomass.data$autumn = ifelse(biomass.data$Treatment %in% c("control","spring"), "no","warming")
#Transform biomass to percentage anomaly
########################################
#get means per species
biomass.data1 = biomass.data %>%
group_by(Species,organ) %>%
summarize(Mean = mean(biomass, na.rm=TRUE))
#merge data
biomass.data1 <-merge(biomass.data, biomass.data1, all.x=T, by=c("Species","organ"))
#get precentage change in biomass
biomass.data1$relativeBiomass = (biomass.data1$biomass / biomass.data1$Mean -1) * 100
#Multivariate linear model
##########################
resultsLM.biomass = biomass.data1 %>%
group_by(Species,organ) %>%
do({model = lm(relativeBiomass ~ spring*autumn, data=.) # create your model
data.frame(tidy(model))}) %>% # get coefficient info
filter(!term %in% c("(Intercept)")) %>% # delete rows
dplyr::select(Species, organ, term, p.value, estimate, std.error) %>% #delete columns
mutate_if(is.numeric, round, 2) %>%
mutate(sig = ifelse(p.value<0.05, "**", ifelse(p.value<=0.1, "*", ""))) #add significance## `mutate_if()` ignored the following grouping variables:
## Columns `Species`, `organ`data.frame(resultsLM.biomass)## Species organ term p.value estimate std.error sig
## 1 Quercus Total springwarming 0.06 35.83 18.03 *
## 2 Quercus Total autumnwarming 0.02 -40.96 17.05 **
## 3 Quercus Total springwarming:autumnwarming 0.55 14.94 24.50
## 4 Quercus Shoot springwarming 0.32 17.43 17.11
## 5 Quercus Shoot autumnwarming 0.01 -43.38 16.18 **
## 6 Quercus Shoot springwarming:autumnwarming 0.20 30.11 23.25
## 7 Quercus Root springwarming 0.02 51.71 21.27 **
## 8 Quercus Root autumnwarming 0.06 -38.87 20.12 *
## 9 Quercus Root springwarming:autumnwarming 0.95 1.86 28.91
## 10 Quercus RSR springwarming 0.06 27.10 13.84 *
## 11 Quercus RSR autumnwarming 0.72 4.81 13.08
## 12 Quercus RSR springwarming:autumnwarming 0.23 -22.94 18.80
## 13 Fagus Total springwarming 0.09 20.95 12.17 *
## 14 Fagus Total autumnwarming 0.38 10.92 12.17
## 15 Fagus Total springwarming:autumnwarming 0.15 -25.67 17.43
## 16 Fagus Shoot springwarming 0.43 9.90 12.51
## 17 Fagus Shoot autumnwarming 0.52 8.06 12.51
## 18 Fagus Shoot springwarming:autumnwarming 0.44 -13.94 17.92
## 19 Fagus Root springwarming 0.02 35.00 14.10 **
## 20 Fagus Root autumnwarming 0.31 14.56 14.10
## 21 Fagus Root springwarming:autumnwarming 0.05 -40.59 20.19 *
## 22 Fagus RSR springwarming 0.03 22.66 9.90 **
## 23 Fagus RSR autumnwarming 0.49 6.94 9.90
## 24 Fagus RSR springwarming:autumnwarming 0.12 -22.38 14.19
## 25 Lonicera Total springwarming 0.81 -1.52 6.14
## 26 Lonicera Total autumnwarming 0.60 3.46 6.49
## 27 Lonicera Total springwarming:autumnwarming 0.98 0.25 8.94
## 28 Lonicera Shoot springwarming 0.35 -6.47 6.77
## 29 Lonicera Shoot autumnwarming 0.62 -3.57 7.16
## 30 Lonicera Shoot springwarming:autumnwarming 0.54 6.05 9.86
## 31 Lonicera Root springwarming 0.42 7.26 8.83
## 32 Lonicera Root autumnwarming 0.10 15.94 9.33 *
## 33 Lonicera Root springwarming:autumnwarming 0.44 -10.06 12.85
## 34 Lonicera RSR springwarming 0.14 14.36 9.37
## 35 Lonicera RSR autumnwarming 0.04 21.13 9.91 **
## 36 Lonicera RSR springwarming:autumnwarming 0.19 -18.31 13.64########
# T-test
########
#get means of control group
###########################
biomass.data1 = biomass.data %>%
filter(Treatment %in% c("control")) %>%
group_by(Species,organ) %>%
summarize(Mean = mean(biomass, na.rm=TRUE))
#merge data
biomass.data1 <-merge(biomass.data[!biomass.data$Treatment=="control",], biomass.data1, all.x=T, by=c("Species","organ"))
#get precentage change in biomass
biomass.data1$relativeBiomass = ifelse(biomass.data1$organ!="RSR", (biomass.data1$biomass / biomass.data1$Mean -1) * 100,
(biomass.data1$biomass - biomass.data1$Mean)*100)
#Ordering of factors
Treatments = c("spring", "autumn", "autumn-spring") #order treatments
biomass.data1$Treatment = factor(biomass.data1$Treatment, levels=Treatments, ordered=T)
Species = c("Quercus", "Fagus", "Lonicera") #order treatments
biomass.data1$Species = factor(biomass.data1$Species, levels=Species, ordered=T)
#T-test
#######
resultsTTbiomass = biomass.data1 %>%
group_by(Species,organ,Treatment) %>%
do({model = t.test(.$relativeBiomass) # create your model
data.frame(tidy(model),
tidy(shapiro.test(x = .$relativeBiomass))[1,2])}) %>%
select(Species, organ, Treatment, p.value, estimate, p.value.1) %>% #delete columns
rename(shapiroTest=p.value.1) %>%
mutate_if(is.numeric, round, 2) %>%
mutate(sig = ifelse(p.value<0.05, "**", ifelse(p.value<=0.1, "*", ""))) #add significance## `mutate_if()` ignored the following grouping variables:
## Columns `Species`, `organ`, `Treatment`#order statistics table to match dataframe
resultsTTbiomass = resultsTTbiomass[with(resultsTTbiomass,
order(resultsTTbiomass$Species, resultsTTbiomass$organ)), ]
data.frame(resultsTTbiomass)## Species organ Treatment p.value estimate shapiroTest sig
## 1 Quercus Total spring 0.08 35.60 0.59 *
## 2 Quercus Total autumn 0.00 -40.68 0.07 **
## 3 Quercus Total autumn-spring 0.40 9.75 0.48
## 4 Quercus Shoot spring 0.25 16.31 0.25
## 5 Quercus Shoot autumn 0.01 -40.60 0.05 **
## 6 Quercus Shoot autumn-spring 0.68 3.89 0.32
## 7 Quercus Root spring 0.06 54.23 0.82 *
## 8 Quercus Root autumn 0.00 -40.76 0.04 **
## 9 Quercus Root autumn-spring 0.33 15.42 0.83
## 10 Quercus RSR spring 0.08 31.47 0.79 *
## 11 Quercus RSR autumn 0.54 5.58 0.59
## 12 Quercus RSR autumn-spring 0.43 10.41 0.65
## 13 Fagus Total spring 0.04 23.08 0.17 **
## 14 Fagus Total autumn 0.08 12.03 0.89 *
## 15 Fagus Total autumn-spring 0.65 6.83 0.07
## 16 Fagus Shoot spring 0.18 10.46 0.56
## 17 Fagus Shoot autumn 0.19 8.51 0.98
## 18 Fagus Shoot autumn-spring 0.79 4.25 0.04
## 19 Fagus Root spring 0.02 40.81 0.21 **
## 20 Fagus Root autumn 0.08 16.97 0.16 *
## 21 Fagus Root autumn-spring 0.48 10.46 0.06
## 22 Fagus RSR spring 0.03 17.95 0.85 **
## 23 Fagus RSR autumn 0.36 5.50 0.97
## 24 Fagus RSR autumn-spring 0.15 5.72 0.02
## 25 Lonicera Total spring 0.70 -1.54 0.55
## 26 Lonicera Total autumn 0.48 3.49 0.67
## 27 Lonicera Total autumn-spring 0.66 2.20 0.14
## 28 Lonicera Shoot spring 0.14 -6.24 0.12
## 29 Lonicera Shoot autumn 0.56 -3.45 0.23
## 30 Lonicera Shoot autumn-spring 0.47 -3.85 0.71
## 31 Lonicera Root spring 0.31 7.96 0.52
## 32 Lonicera Root autumn 0.07 17.49 0.38 *
## 33 Lonicera Root autumn-spring 0.05 14.41 0.78 *
## 34 Lonicera RSR spring 0.09 8.20 0.23 *
## 35 Lonicera RSR autumn 0.05 12.06 0.76 *
## 36 Lonicera RSR autumn-spring 0.01 9.81 0.72 **5. Figures
Figure 1
#2D density plot
Dens2D = ggplot(clim.data, aes(x=outside_tmp, y=inside_tmp-outside_tmp) ) +
stat_density_2d(aes(fill = ..level..), geom = "polygon") +
scale_fill_distiller(palette=4, direction=1) +
geom_hline(aes(yintercept=4.7))+
coord_cartesian(ylim = c(0, 10),xlim = c(-4, 19))+
xlab("Ambient temperature (°C)") +
ylab("Deviation from control (°C)") +
plotTheme
#Boxplot
boxp = ggplot(clim.data2, aes(x=treatment, y=diffOutIn, fill=treatment)) +
geom_boxplot(outlier.shape=NA,alpha=.6) +
coord_flip(ylim = c(0.5, 10)) +
xlab("") +
ylab("") +
scale_fill_manual(values=c("orange","green4"))+
plotTheme2
#Histogram
histo = ggplot(clim.data2, aes(x=diffOutIn, fill=treatment)) +
geom_histogram(binwidth=.5, alpha=.6, position="identity") +
geom_vline(data=clim.data1, aes(xintercept=medianDiff, colour=c("green4","orange3"),alpha=.8),
linetype="dashed", size=1)+
xlab("Deviation from control (°C)") +
ylab("Count (number of days)") +
scale_fill_manual(values=c("orange","green4"))+
scale_colour_manual(values=c("orange","green3"))+
coord_cartesian(xlim = c(0.5, 10), ylim = c(3, 60))+
plotTheme3
#Time plot
#extract date information
phenoMean$dateMin = as.Date(phenoMean$Mean-phenoMean$SD, origin = paste0(phenoMean$Year,"-01-01"))
phenoMean$dateMax = as.Date(phenoMean$Mean+phenoMean$SD, origin = paste0(phenoMean$Year,"-01-01"))
phenoMean$Y = ifelse(phenoMean$Species=="Lonicera",30,ifelse(phenoMean$Species=="Quercus",28,26))
TempPlot = ggplot() +
geom_rect(data=phenoMean[phenoMean$phenology=='Leaf_out',],
mapping=aes(xmin=dateMin, xmax=dateMax, ymin=-Inf, ymax=Y),
fill=rep(c('green4','green1',"green3"),each=3), alpha=0.4) +
geom_rect(data=phenoMean[phenoMean$phenology=='Senescence',],
mapping=aes(xmin=dateMin, xmax=dateMax, ymin=-Inf, ymax=Y),
fill=rep(c('orange4','orange1',"orange3"),each=3), alpha=0.4) +
geom_line(data=clim.data3, aes(x = date, y= value, group = variable, colour = variable2)) +
geom_segment(data=phenoMean[phenoMean$phenology=='Leaf_out',],
aes(x = dateMin, y = Y, xend = dateMax, yend = Y),
size=3,
col=rep(c('green4','green1',"green3"),each=3) )+
geom_segment(data=phenoMean[phenoMean$phenology=='Senescence',],
aes(x = dateMin, y = Y, xend = dateMax, yend = Y),
size=3,
col=rep(c('orange4','orange1',"orange3"),each=3) )+
scale_color_manual(values=c("orange3","black","green4"))+
coord_cartesian(ylim=c(-7,29))+
xlab("") +
ylab("Temperature (°C)") +
plotTheme
#define plot layout
layout <- "
AA
AA
AA
BD
CD
CD"
#Merge plots
Fig1_ClimPlot = TempPlot + boxp + histo + Dens2D + plot_layout(design = layout)
Fig1_ClimPlot
##########
#Save PDFs
##########
pdf(paste(output.dir,"Fig1_Temperature deviation.pdf",sep="/"), width=9, height=6, useDingbats=FALSE)
Fig1_ClimPlot
dev.off()## quartz_off_screen
## 2Figure 2
#create panel labels
dat_text = data.frame(
label = c("a","b","c","d","e","f","g","h","i"),
phenology = rep(c("Leaf_out","Senescence","GSL"),3),
Species = c(rep("Quercus",3), rep("Fagus",3), rep("Lonicera",3)),
term = "springwarming")
#Plot
Fig2_LMplot = ggplot(data = resultsLM.pheno, mapping = aes(x = term, y = estimate,
fill=phenology, alpha=term)) +
geom_bar(position=position_dodge(), stat="identity", color="black") +
geom_errorbar(aes(ymin=estimate-std.error, ymax=estimate+std.error),
width=.2, # Width of the error bars
position=position_dodge(.9))+
geom_hline(yintercept=0) +
geom_text(data = dat_text, mapping = aes(x = Inf, y = Inf,
hjust = 1, vjust = 1,
label = label,
fontface = "bold"))+
stat_summary(geom = 'text', label = resultsLM.pheno$sig,
fun.y = mean, vjust = 2) +
scale_alpha_discrete(range = c(1, 0.7)) +
xlab("Phenophase") +
ylab("Effect size (days)") +
scale_fill_manual(values = c("green3", "orange", "blue")) +
facet_grid(Species~phenology) +
plotTheme4
Fig2_LMplot
#########
#Save PDF
#########
pdf(paste(output.dir,"Fig2_PhenologyLinearTreatmentModel.pdf",sep="/"), width=7, height=7, useDingbats=FALSE)
Fig2_LMplot
dev.off()Figure 3
volume.data$Species = factor(volume.data$Species, levels=c("Quercus robur", "Fagus sylvatica", "Lonicera xylosteum"))
pd = position_dodge(0.2)
Fig3_LinePlot = ggplot(volume.data, aes(x=Year, y=Volume, colour=Treatment, group=Treatment)) +
stat_summary(fun = mean, geom="line", position=pd) +
stat_summary(fun.data = "mean_se", size = 0.5, position=pd) +
labs(x = "Year", y = "Stem volume (cm3)") +
scale_color_manual(values=c("green3", "orange", "red3", "black")) +
coord_cartesian(ylim = c(8, 152))+
facet_grid(~Species) +
plotTheme
Fig3_LinePlot
#########
#Save PDF
#########
pdf(paste(output.dir,"Fig3_LinePlot.pdf",sep="/"), width=9, height=4, useDingbats=FALSE)
Fig3_LinePlot
dev.off()Figure 4
#Ordering of factors
term = c("springwarming", "autumnwarming", "springwarming:autumnwarming")
resultsLM.biomass$term = factor(resultsLM.biomass$term, levels=term, ordered=T)
#create panel labels
dat_text = data.frame(
label = c("a","b","c","d","e","f","g","h","i","j","k","l"),
organ = rep(c("Total","Shoot","Root","RSR"),3),
Species = c(rep("Quercus",4), rep("Fagus",4), rep("Lonicera",4)),
term = "springwarming")
#Plot
LMplot = ggplot(data = resultsLM.biomass, mapping = aes(x = term, y = estimate,
fill=organ, alpha=term)) +
geom_bar(position=position_dodge(), stat="identity", color="black") +
geom_errorbar(aes(ymin=estimate-std.error, ymax=estimate+std.error),
width=.2, # Width of the error bars
position=position_dodge(.9))+
geom_hline(yintercept=0) +
geom_text(data = dat_text, mapping = aes(x = -Inf, y = -Inf,
hjust = -0.3, vjust = -1,
label = label,
fontface = "bold"))+
stat_summary(geom = 'text', label = resultsLM.biomass$sig,
fun.y = mean, vjust = 2) +
scale_fill_manual(values=c("orange","yellow2","red3","grey70"))+
scale_alpha_discrete(range = c(1, 0.7)) +
xlab("Phenophase") +
ylab("Effect size (%)") +
facet_grid(Species~organ) +
plotTheme4
LMplot
#########
#Save PDF
#########
pdf(paste(output.dir,"Fig4_BiomassLinearTreatmentModel.pdf",sep="/"), width=8, height=7, useDingbats=FALSE)
LMplot
dev.off()Figure 5
#Effect sizes
#############
#create panel labels
dat_text = data.frame(
label = c("a","b","c","d","e","f","g","h","i","j","k","l"),
organ = rep(c("Total","Shoot","Root","RSR"),3),
Species = c(rep("Quercus",4), rep("Fagus",4), rep("Lonicera",4)),
phenology = "Leaf_out")
#Plot
LMplot = ggplot(data = resultsLM.pheno2, mapping = aes(x = phenology, y = estimate,
fill=organ, alpha=phenology)) +
geom_bar(position=position_dodge(), stat="identity", color="black") +
geom_errorbar(aes(ymin=estimate-std.error, ymax=estimate+std.error),
width=.2, # Width of the error bars
position=position_dodge(.9))+
geom_hline(yintercept=0) +
geom_text(data = dat_text, mapping = aes(x = -Inf, y = -Inf,
hjust = -0.1, vjust = -1,
label = label,
fontface = "bold"))+
stat_summary(geom = 'text', label = resultsLM.pheno2$sig,
fun.y = mean, vjust = 2) +
scale_fill_manual(values=c("orange","yellow2","red3","grey70"))+
scale_alpha_discrete(range = c(1, 0.7)) +
xlab("Phenophase") +
ylab("Effect size (%/day)") +
facet_grid(Species~organ) +
plotTheme4
LMplot
#Biomass
########
#order
Organs = c("Total","Shoot","Root","RSR") #order treatments
pheno.data3$organ = factor(pheno.data3$organ, levels=Organs, ordered=T)
#create panel labels
dat_text = data.frame(
label = c("a","b","c","d","e","f","g","h","i"),
phenology = rep(c("Leaf_out","Senescence","GSL"),3),
Species = c(rep("Quercus",3), rep("Fagus",3), rep("Lonicera",3)),
organ = "Root")
BiomassPlot = ggplot(pheno.data3[pheno.data3$organ %in% c("Total","Shoot","Root"),],
aes(x=phenologyAnomaly, y=relativeBiomass, colour=organ)) +
geom_hline(yintercept=0) +
geom_vline(xintercept=0) +
geom_smooth(method=lm, fullrange = T, alpha = 0.2) +
labs(x = "Phenology anomaly (days)", y = "Biomass anomaly (%)") +
scale_color_manual(values=c("orange","yellow2","red3"))+
geom_text(data = dat_text, mapping = aes(x = -Inf, y = Inf,
hjust = -0.2, vjust = 1.3,
label = label,
fontface = "bold"), color="black")+
geom_text(data = resultsLM.pheno2[resultsLM.pheno2$organ=="Total",],
mapping = aes(x = Inf, y = Inf, hjust = 1.5, vjust = 2.5,
label = paste(round(estimate,1), " %/day", sig, sep="")),
size=3.5, color="orange")+
geom_text(data = resultsLM.pheno2[resultsLM.pheno2$organ=="Shoot",],
mapping = aes(x = Inf, y = Inf, hjust = 1.5, vjust = 4.5,
label = paste(round(estimate,1), " %/day", sig, sep="")),
size=3.5, color="yellow2")+
geom_text(data = resultsLM.pheno2[resultsLM.pheno2$organ=="Root",],
mapping = aes(x = Inf, y = Inf, hjust = 1.5, vjust = 6.5,
label = paste(round(estimate,1), " %/day", sig, sep="")),
size=3.5, color="red3")+
facet_grid(Species~phenology, scales = "free") +
plotTheme
BiomassPlot
#RSR
RSRplot = ggplot(pheno.data3[pheno.data3$organ %in% c("RSR"),], aes(x=phenologyAnomaly, y=biomass)) +
geom_point(size=0.2) +
geom_smooth(method=lm, fullrange = F) +
labs(x = "Day", y = "Root:shoot ratio") +
scale_color_viridis(option="viridis",discrete=TRUE) +
geom_text(data = resultsLM.pheno2[resultsLM.pheno2$organ=="RSR",],
mapping = aes(x = Inf, y = Inf, hjust = 1.5, vjust = 2.5,
label = paste("Slope = ", round(estimate,3), " day-1", sig, sep="")),
size=3.5, color="black")+
facet_grid(Species~phenology, scales = "free") +
plotTheme
RSRplot
#########
#Save PDF
#########
pdf(paste(output.dir,"Fig5_PhenologyBiomass.pdf",sep="/"), width=8, height=7, useDingbats=FALSE)
BiomassPlot
RSRplot
dev.off()Figure S1
#create panel labels
dat_text <- data.frame(
label = c("a","b","c","d","e","f","g","h","i"),
phenology = rep(c("Leaf_out","Senescence","GSL"),3),
Species = c(rep("Quercus",3), rep("Fagus",3), rep("Lonicera",3)),
Treatment = "spring")
#Plot
FigS1_PhenologyChangePlot = ggplot(data = pheno.data2, mapping = aes(x = Treatment, y = phenologyChange,
fill=phenology, alpha=Treatment)) +
stat_summary(fun.y = mean,
geom = "bar", color="black") +
stat_summary(fun.data = mean_se,
geom = "errorbar", width=0.2) +
geom_hline(yintercept=0) +
geom_text(data = dat_text, mapping = aes(x = -Inf, y = Inf,
hjust = -0.2, vjust = 1.5,
label = label,
fontface = "bold"))+
stat_summary(geom = 'text', label = resultsTT$significance,
fun.y = mean, vjust = 2) +
scale_alpha_discrete(range = c(1, 0.7)) +
xlab("Treatment") +
ylab("Phenological change (days)") +
scale_fill_manual(values = c("green3", "orange", "blue")) +
facet_grid(Species~phenology) +
plotTheme4
FigS1_PhenologyChangePlot
#########
#Save PDF
#########
pdf(paste(output.dir,"FigS1_PhenologyChange.pdf",sep="/"), width=7, height=7, useDingbats=FALSE)
FigS1_PhenologyChangePlot
dev.off()Figure S2
#Plot
#####
#create panel labels
dat_text <- data.frame(
label = c("a","b","c","d","e","f","g","h","i","j","k","l"),
organ = rep(c("Total","Shoot","Root","RSR"),3),
Species = c(rep("Quercus",4), rep("Fagus",4), rep("Lonicera",4)),
Treatment = "spring")
#Plot
BiomassChangePlot = ggplot(data = biomass.data1, mapping = aes(x = Treatment, y = relativeBiomass,
fill=organ, alpha=Treatment)) +
stat_summary(fun.y = mean,
geom = "bar", color="black") +
stat_summary(fun.data = mean_se,
geom = "errorbar", width=0.2) +
geom_hline(yintercept=0) +
geom_text(data = dat_text, mapping = aes(x = -Inf, y = -Inf,
hjust = -0.3, vjust = -1,
label = label,
fontface = "bold"))+
stat_summary(geom = 'text', label = resultsTTbiomass$sig,
fun.y = max, vjust = 2.9) +
scale_alpha_discrete(range = c(1, 0.7)) +
coord_cartesian(ylim = c(-60, 60))+
xlab("Treatment") +
ylab("Biomass change (%)") +
scale_fill_manual(values=c("orange","yellow2","red3","grey70"))+
facet_grid(Species~organ) +
plotTheme4
BiomassChangePlot
#########
#Save PDF
#########
pdf(paste(output.dir,"FigS2_BiomassChange.pdf",sep="/"), width=8, height=7, useDingbats=FALSE)
BiomassChangePlot
dev.off()6. Reproducibility
Reproducibility receipt
## datetime
Sys.time()## [1] "2021-04-28 11:10:32 CEST"## repository
#git2r::repository()
## sessioninfo
sessionInfo()## R version 3.6.2 (2019-12-12)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] car_3.0-6 carData_3.0-3 viridis_0.5.1 viridisLite_0.3.0
## [5] broom_0.5.4 stringr_1.4.0 git2r_0.27.1 data.table_1.12.8
## [9] lubridate_1.7.4 patchwork_1.0.0 dplyr_0.8.4 ggplot2_3.2.1
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.0.0 xfun_0.12 reshape2_1.4.3 purrr_0.3.3
## [5] haven_2.2.0 lattice_0.20-38 colorspace_1.4-1 generics_0.0.2
## [9] vctrs_0.2.2 htmltools_0.4.0 yaml_2.2.1 rlang_0.4.4
## [13] pillar_1.4.3 foreign_0.8-72 glue_1.3.1 withr_2.1.2
## [17] RColorBrewer_1.1-2 readxl_1.3.1 plyr_1.8.5 lifecycle_0.1.0
## [21] cellranger_1.1.0 munsell_0.5.0 gtable_0.3.0 zip_2.0.4
## [25] evaluate_0.14 labeling_0.3 knitr_1.28 rio_0.5.16
## [29] forcats_0.4.0 curl_4.3 Rcpp_1.0.3 scales_1.1.0
## [33] backports_1.1.5 abind_1.4-5 farver_2.0.3 gridExtra_2.3
## [37] hms_0.5.3 digest_0.6.23 stringi_1.4.5 openxlsx_4.1.4
## [41] grid_3.6.2 tools_3.6.2 magrittr_1.5 lazyeval_0.2.2
## [45] tibble_2.1.3 crayon_1.3.4 tidyr_1.0.2 pkgconfig_2.0.3
## [49] MASS_7.3-53 assertthat_0.2.1 rmarkdown_2.1 R6_2.4.1
## [53] nlme_3.1-142 compiler_3.6.2